gusucode.com > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM源码程序 > 支持向量机工具箱 - LIBSVM OSU_SVM LS_SVM\stprtool\svm\kerskfmat.m
function [Alpha,bias,sol,t,flps,margin,up,lo]=... kerskfmat(X,I,epsilon,ker,arg,tmax,C) % KERSKFMAT fast version of KERNELSK. % [Alpha,bias,sol,t,flps,margin,up,lo]=kerskfmat(X,I,epsilon,ker,arg,tmax,C) % % Faster version (with respect to number of floating point operations) % of KERNELSK algorithm. It does not compute the whole kernel matrix. % % Inputs: % X [NxL] training patterns, N is dimension and L number of patterns. % I [1xL] labels, 1 for 1st class and 2 for 2nd class. % epsilon [1x1] precision of found solution. The margin of found % hyperplane is less than the optimal margin at most by epsilon. % ker [string] kernel, see 'help kernel'. % arg [...] argument of given kernel, see 'help kernel'. % tmax [int] maximal number of iterations. % C [real] trade-off between margin and training error. % % Outputs: % Alpha [1xL] Weights (Lagrangians) of patterns. % bias [real] bias (threshold) of found decision rule. % sol [int] 1 solution is found % 0 algorithm stoped (t == tmax) before converged. % -1 hyperplane with margin greater then epsilon % does not exist. % t [int] number of iterations. % margin [real] margin between classes. % flps [int] number of used floating point operations. % up [1,t] evolution of the upper bound on the optimal margin. % lo [1,t] evolution of the lower bound on the optimal margin. % CE [real] classification error on training patterns. % % See also KERNELSK, SVM. % % Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac % (c) Czech Technical University Prague, http://cmp.felk.cvut.cz % Written Vojtech Franc (diploma thesis) 02.11.1999, 13.4.2000 % Modifications % 19-September-2001, V.Franc, comments changed. % 18-August-2001, V.Franc, up and lo bounds added % 9-August-2001, V.Franc, version without computing kernel matrices % 13-July-2001, V.Franc, comments % 12-July-2001, V.Franc, C, bias and normal vect. normalized. % 11-July-2001, V.Franc, Rosta Horcik proved that the computation % of threshold is OK. % 10-July-2001, V.Franc, derived from kekozinec2 flops(0); % set default values of the input argiments if nargin < 7, C = inf; end % maximal number of iteraions if nargin < 6, tmax = inf; end % indexes of pattens in the 1st and 2nd class xinx1 = find(I == 1); xinx2 = find(I == 2); X1=X(:,xinx1); % patters from 1st class X2=X(:,xinx2); % patters from 1st class l1 = size(X1,2); % number os patterns l2 = size(X2,2); % make problem lin-separable in high dimensional space if C ~= 0, kadd = 1/(2*C); else kadd = 0; end % convex coeficients defining normal of the decision hyperplane % (they correspond to the Lagrangian multiplyers). s1 = zeros(l1, 1); s2 = zeros(l2, 1); % initial sol s1(1)=1; % take the 1st pattern from the 1st class s2(1)=1; % take the 2nd pattern from the 2st class sol=0; t = 0; minXDA1=-inf; minXDA2=-inf; % -- Inicalization of temp. variables --------------------------------- Di1 = zeros(l1,1); Fi1 = zeros(l1,1); for i=1:l1, Di1(i) = kernel(X1(:,1),X1(:,i), ker,arg); Fi1(i) = kernel(X2(:,1),X1(:,i), ker,arg); end Di1(1) = Di1(1) + kadd; Di2 = zeros(l2,1); Fi2 = zeros(l2,1); for i=1:l2, Di2(i) = kernel(X1(:,1),X2(:,i), ker,arg); Fi2(i) = kernel(X2(:,1),X2(:,i), ker,arg); end Fi2(1) = Fi2(1)+kadd; a = kernel(X1(:,1),X1(:,1),ker,arg) + kadd; b = kernel(X2(:,1),X2(:,1),ker,arg) + kadd; c = kernel(X1(:,1),X2(:,1),ker,arg ); % upper and lower bounds on the optimal margin up=[]; lo=[]; % main cycle while sol == 0 & tmax > t, t = t + 1; sol = 1; % -- compute auxciliary variables -- if sqrt( a -2*c +b) <= 0, % algorithm has converged to the zero vector --> classes overlap sol = -1; break; end [emin,inx1] = min(Di1-Fi1); [gmin,inx2] = min(Fi2-Di2); % projection x \in X_1 on (w_1 - w_2) proj1 = (emin + b -c )/sqrt(a-2*c+b); % projection x \in X_2 on (w_2 - w_1) proj2 = (gmin + a - c)/sqrt(a-2*c+b); % --- compute stop condition for the alpha1 (1st class) ------ % (proj1 < proj2) ~ the worst point will be used for update if (proj1 < proj2) & (proj1 <= (sqrt(a-2*c+b) - epsilon)), % -- Adaptation phase of vector alpha1 ---------------------------- k = (a - emin - c)/... (a+kernel(X1(:,inx1),X1(:,inx1),ker,arg)+kadd-2*(Di1(inx1)-Fi1(inx1)) ); k = min( 1, k ); s1 = s1 * (1-k); s1(inx1) = s1(inx1) + k; sol = 0; % ------------------------------------------------------------- a = a*(1-k)^2 + 2*(1-k)*k*Di1(inx1) + ... k^2 * (kernel(X1(:,inx1),X1(:,inx1),ker,arg)+kadd ); c = c*(1-k) + k*Fi1(inx1); for i=1:l1, Di1(i) = Di1(i)*(1-k) + k*kernel(X1(:,i),X1(:,inx1),ker,arg); end Di1(inx1) = Di1(inx1) + k*kadd; for i=1:l2, Di2(i) = Di2(i)*(1-k) + k*kernel(X2(:,i),X1(:,inx1),ker,arg); end else % --- compute stop condition for the alpha2 (2st class) ------ if proj2 <= (sqrt(a-2*c+b) - epsilon ), % -- Adaptation phase ---------------------------------- k = (b - gmin -c)/... (b+kernel(X2(:,inx2),X2(:,inx2),ker,arg)+kadd-2*(Fi2(inx2)-Di2(inx2))); k = min( 1, k ); s2 = s2 * (1-k); s2(inx2) = s2(inx2) + k; sol = 0; % ------------------------------------------------------ b = b*(1-k)^2 + 2*(1-k)*k*Fi2(inx2) + ... k^2 * (kernel(X2(:,inx2),X2(:,inx2),ker,arg)+kadd ); c = c*(1-k) + k*Di2(inx2); for i=1:l1, Fi1(i) = (1-k)*Fi1(i) + k*kernel(X2(:,inx2),X1(:,i),ker,arg); end for i=1:l2, Fi2(i) = (1-k)*Fi2(i) + k*kernel(X2(:,inx2),X2(:,i),ker,arg); end Fi2(inx2) = Fi2(inx2) + k*kadd; end end % store up=||w||/2 and current margin m(w1-w2,theta) = min( m1, m2) m = min([proj1,proj2]) - 0.5*sqrt(a-2*c+b); up = [up,sqrt(a-2*c+b)/2]; lo = [lo,m ]; %% disp(sprintf('step = %d', t )); end if sol == 1 & (proj1 < 0 | proj2 < 0), sol = -2; end % -- compute threshold ----------------------- % sqared margin in transfromed space margin2 = a - 2*c + b; % threshold after normalization theta = (a-b)/margin2; % solution (normal vect. in the transformed space) after normalization s1=2*s1/margin2; s2=2*s2/margin2; % -- make SVM-like output -------------------- Alpha=zeros(1,l1+l2); Alpha(xinx1)=s1; Alpha(xinx2)=s2; bias = -theta; % -- compute margin ----------------------------------- if nargout >= 6, margin = 0; for i=find(Alpha ~= 0), for j=find(Alpha ~= 0 ), margin = margin + Alpha(i)*Alpha(j)*itosgn(I(i))*... itosgn(I(j))* kernel(X(:,i),X(:,j),ker,arg); end end margin = 2/sqrt(margin); end % overall number of used float point operations flps=flops; return;